function kbasis = make_kbasis_2(varargin)
%
% MAKE_KBASIS makes kbasis matrix or structure for efficient computation
%
% MAKE_KBASIS(FID,MODEL,METHOD[,NCPUS])

% MAKE_KBASIS always expects:
%    * FID.KX -- kx-coordinate vector
%    * FID.KY -- ky-coordinate vector
%    * FID.XX -- actual x (and y) coordinates in cm
%
% If METHOD='separable', MAKE_KBASIS expects:
%    * FID.IX --  vector of ROS row-coordinate indices
%    * FID.IY --  vector of ROS col-coordinate indices
%    * FID.RESO -- image side size in voxels
% and returns the following structure:
%    * KBASIS.EKXX -- exp(j2pi kx xx) matrix
%    * KBASIS.EKYY -- exp(j2pi ky yy) matrix
%    * KBASIS.IX = FID.IX
%    * KBASIS.IY = FID.IY
%
% If METHOD~='separable', MAKE_KBASIS expects:
%    * FID.X --  vector of ROS x-coordinates
%    * FID.Y --  vector of ROS y-coordinates
% Optional:
%    * MODEL.N -- number of voxels in ROS
%    * FID.N -- number of FID samples
% and returns a matrix of gradient responses (col) for each (x,y) (row).
%
% NCPUS indicates the number of CPU's for parallel implementation

coex=-2i*pi;

fid = varargin{1};
model = varargin{2};
method = varargin{3};
if nargin==4,
    Ncpus = varargin{4};
else
    Ncpus = 1;
end

if isfield(fid,'N') 
    NF = fid.N;
else
    NF = length(fid.sig1);
end
if isfield(model,'N') 
    NM = model.N; 
else
    NM = length(model.pvec)/2; 
end
% precalculate arrays to make fid synthesis more efficient
if (strcmp(method,'separable'))
    % separable approach
    if (Ncpus < 2) % serial computing setup
        kbasis.ekxx = exp(fid.kx*(coex*model.ix1d));
        kbasis.ekyy = exp(fid.ky*(coex*model.ix1d));
    else % parallel computing setup
        for ii = 1:Ncpus
            kbasis.ekxx(:,:,ii) = exp(fid.kx(ii,:)'*(coex*model.ix1d));
            kbasis.ekyy(:,:,ii) = exp(fid.ky(ii,:)'*(coex*model.ix1d));
        end
    end
    kbasis.iy = model.iy;
    kbasis.ix = model.ix;
    kbasis.ix1d = 1:model.reso;
else
    % kbasis is the array of ideal local responses to the encoding gradients
    % for unit-magnitude on-resonance magnetization vectors at each pixel
    % (x,y) location.  These are "basis" functions of a sort.
    if (Ncpus < 2) % serial computing setup

%         for nl = 1:NM
%             kbasis(nl,:) = exp( coex * ( tempkx * model.x(nl) + tempky * model.y(nl) ) );
%         end
    clear kbasis
    kbasis = single(exp( coex * ( model.x(1:NM) * fid.kx(1:NF) + model.y(1:NM) * fid.ky(1:NF))));
    else % parallel computing setup
        temp = zeros(NM,NF);
        kbasis = zeros(NM,NF/Ncpus,Ncpus);
        for nl = 1:NM
            temp(nl,:) = exp(coex*(fid.kx(1:NF)*model.x(nl)+fid.ky(1:NF)*model.y(nl)));
        end
        for ii = 1:Ncpus
            kbasis(:,:,ii) = temp(:,ii:Ncpus:NF);
        end
    end
end


